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Abstract 

We have proposed a cluster heat bath method in Monte Carlo simulations 

of Ising models in which one of the possible spin configurations of a cluster is 

selected in accordance with its Boltzmann weight. We have argued that the 

method improves slow relaxation in complex systems and demonstrated it in 

an axial next-nearest-neighbor Ising(ANNNI) model in two-dimensions. 
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Recently various algorithms have been proposed to reduce the CPU time of computer in 
the Monte Carlo(MC) simulation [|I|. Cluster-flip algorithms 0,^ were proposed using ideas 
from percolation theory [^]. Although the methods were very efficient to simulate large sys- 
tems near criticality, these were not successfully applied to complex systems which contain 
frustrated interactions such as spin-glasses. On the other hand, to study the ordered state of 
complex systems, extended ensemble methods were developed . In these methods, how- 
ever, the conventional single-spin-fiip algorithm is used to guarantee the ergodicity. Quite 
recently, a new update method was proposed in which the spin configuration of a chain of 
Ising spins is updated in accordance with the Boltzmann weight [H- The method is very 



effective for quasi-one-dimensional models and enables us to make realistic simulations |10 



of quasi-one-dimensional Ising magnets such as CsCoBra [|ri|] and CsCoCla [|T2|. However, 
it was not so useful for ordinary two- and three-dimensional models, especially for complex 
systems. 

In this Letter, we propose a general configuration-update algorithm which is applicable 
to various Ising models with short-range interactions and very effective for improving slow 
relaxation in complex systems. We consider to find the probable spin configuration of a 
cluster of N spins. This can be readily done when is small, but becomes difficult when is 
increased. However, if the cluster is decomposed into layers of spins and interlayer couplings 
exist only between the spins on adjacent layers, we can update the spin configuration of 
the layers step by step with the aid of a transfer matrix technique. Thus, we can treat a 
larger cluster, e.g., a cluster oi N ^ M x L spins in a cubic lattice with the nearest neighbor 
couplings, where L is the linear size of the lattice and M is the number of the spins of the 
layer. We can prove that the spin configuration realized in this method is in accordance with 
the Boltzmann weight. So we call the algorithm a cluster heat bath (CHB) method. We 
argue that the CHB method improves slow relaxation in complex systems and demonstrate 
it in an axial next-nearest-neighbor Ising (ANNNI) model in two- dimensions. 

We start with an Ising model described by the Hamiltonian 
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n^- Jij(^i(^j^ (1) 

<i,j> 

where ai{— ±1) are Ising spins and Jij are coupling constants. Now we pick out a cluster of 
spins and consider its probable spin configuration under the condition that the other spins 
are fixed. It is noted that the cluster defined here is an ensemble of connected spins in a 
fixed part of the system, not one of connected spins with the same sign. We suppose that 
the cluster is composed of L layers and interactions exist only between neighboring layers. 
Hereafter we call the cluster as Al- The Hamiltonian of Ax, is, then, described as 

L-l L 

'^Al = - E ^1,1+^ -^Bi, (2) 
1=1 1=1 



with 



»-,Hi=i:i:4'""-f-r'. (3) 



where jf^ and J-j'^^^^ are exchange interactions between ith and jth spins on the same layer 
(l) and those between different layers (l) and (Z + 1), respectively. The effective field hf^ of 
ith spin on the layer (/) is given as a sum of the external field and the exchange field which 
comes from the spins surrounding A^: 

^Y.Ji3^3+^H, (5) 
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where m and H are the magnetic moment and the external field, respectively. 

Now we consider the probable spin configuration {crl^-*} of the layer (L). We define the 
following weight function Fi{{af''}): 

Fi{WT})= E ••• E exp(/?X:i/,,,+i + /?X:i^.), (6) 

where (3 — l/ksT with ks and T being the Boltzmann's constant and temperature, respec- 
tively. This function can be readily obtained from the recursion formula 

HW'P})- E F,_,{{at'^})eMm-i,i + m l>2 (7) 
{^f-i)=±i} 



with the initial function Fi({al^^}) — exp(/3Si). The probabihty Pl({(tI^^}) of the spin 
configuration {crl^^} of the layer (L) is given as 

where Zl{— J2^^iL)__^^y -^-^({'^i^''})) ^^e partition function of Al. Thus, we can determine 
the spin configuration {crl^^} of the layer (L) using a uniform random number. 

The next step is to determine the spin configuration of the layer (L—l) under the 

condition that the spin configuration of the layer (L) is given as {cr^^^}. This is equivalent to 
determine the spin configuration {a^^ of the cluster Al-i which is obtained by removing 
the layer (L) from A^. The spins on this layer (L) of Al now contribute to the effective 
fields on the layer (L — 1) of A^-i- Thus, the effective field hi^~^^ of the ith spin on the 
layer (L — 1) of Al-i is given as 

}^-^^h{t-^+Y.jt''''-f- (9) 
j 

Then, the function FL_i{{a^^~^^}) oi A^ becomes 

FL_i({af = ^ Fi_2({af exp{PHL-2,L-i + PBl-i + I3Hl-i,l) 

= F^_i({af ^^)}) cxp(/3i/z.-i,L) (10) 

for Al-i. The probability PL-i{{af'~^'^}) of the spin configuration {af' ^•'} of the layer 
(L — 1) is given as 

= (11) 



with 



= Fa{aP})exp(-/55^). (12) 

Thus, we can also determine the spin configuration {cr-^"^^} of the layer (L — l). Repeating 
this procedure from the layer (L — 1) to the layer (1), the spin configuration of the cluster 
Al can be updated. 



We can readily show that the spin configuration of Al reahzed in this procedure 
is in accordance with the Boltzmann weight. The probabihty of this spin configura- 
tion PA^{{(Ji^^}, {a },■■■, {al^^}) is given as the product of the individual probabilities 
PiiWh), i.e., 

1=1 

1=1 A 

1 



exp{-f3nAj. (13) 



This is nothing but the Boltzmann weight. Hence we call this algorithm a cluster heat bath 
(CHB) method. When all the spins are updated one time, we call it one MC sweep just like 
in the conventional MC method. 

The CHB method can be applied to clusters with any shape. It is not necessary that 
the interactions of the model are of the nearest neighbor. Only restriction is that we can 
decompose the cluster into layers as described in the form of eq. (2). It is most effective 
to choose the cluster as ladders with its width of M in two-dimensions and columns with 
its intersection of x My in three-dimensions, because we can use the transfer matrix 
technique [|T^. The CHB method may improve relaxation in various systems. The updated 
spin configuration of the cluster depends on its environment but not on its original spin 
configuration. Thus, the spin structure may always fiuctuate in the scale of the cluster size 
and, of course, a cluster-fiip effect is automatically taken into account W^. Moreover, if we 



choose the cluster appropriate to the model, we may perform most effective MC simulation. 



Here, we demonstrate it in a simulation of the ANNNI model |]13| in two- dimensions. 

The ANNNI model is an array of Ising chains with ferromagnetic interaction Ji between 
spins in adjacent chains and competing antiferro magnetic interaction J2 between spins in 
next-nearest-neighbor chains, augmented by ferromagnetic nearest-neighbor interaction Jq 
in the chains. It is well known that, for k = — J2/J1 > \, the ground state of the model is 
the (2, 2) antiphase described by an alternate arrangement of two up-spin chains and two 



down-spin chains, i.e., ■ ■ ■ + H !-+■■■. It is believed that a floating incommensurate (IC) 

phase appears between the (2, 2) antiphase and the paramagnetic phase. However, it turned 
out to be much CPU time comsuming to estabhsh the phase boundaries rehably |13[. The 
difficulty in the MC simulation of the ANNNI model is that the spin structure depends on 
initial spin configurations. This is because the spin structure of the IC phase near transition 

to the (2, 2) antiphase consists of regions of the (2, 2) antiphase separated by + + + or 

walls |T^. Therefore, it is necessary to insert 4 walls simultaneously to the (2,2) antiphase 
to get the IC phase starting from the (2, 2) antiphase. That is, we must rearrange at least 
16Nx spins to get the IC phase, where is the number of the spins of the chain, which 



is not easy to be realized by using the conventional single-spin-ffip MC method [0. This 
difficulty is not largely relieved even when we choose an open boundary condition, because 



open boundaries lead to a pinning effect |T^, i.e., the end two chains tend to take either ++ 



or configuration. In this case, we must rearrange at least SA^^: spins at the ends. So the 



ANNNI model is one of the most difficult models in the computer simulation |T^. However, 
if we use the CHB method with clusters of M x A^. spins with M > 16 (or M > 8 at the 
ends), we may easily add or remove the walls. 

To examine our speculation, we performed the CHB simulation of the model with Jq = Ji 
and K = 0.6 on the lattice of x Ny = 64 x 128 spins with open boundary conditions. 



We treated clusters of 8 x 64 spins |jl9[. At each temperature, we started with two different 
initial spin configurations, i.e., a paramagnetic spin configuration and the (2,2) antiphase 
spin configuration, and calculated quantities of interest. Here we present results of the 
square of the chain magnetization M2 which plays the role of the order parameter of this 



model 20,21 



1 ^'a 1 Afo; 

M-^B^E-.)'- (14) 

As temperature is decreased below T = 1.0 Ji, the relaxation becomes very slow. We present, 
in Figs. 1(a) and 1(b), typical results of the MC sweep dependence of M2 at T = 0.9Ji in 
the conventional MC and CHB methods, respectively. Figure 1(a) shows most clearly the 
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difficulty of the computer simulation of the ANNNI model as mentioned above. However, 
we could get its equilibrium value within a reasonable number of MC sweeps using the CHB 
method as seen in Fig. 1(b). We calculated the average values (M2) at different temperatures 
using both the methods. Results are presented in Figs. 2(a) and 2(b) for the conventional 
MC method and the CHB method, respectively. At all temperatures, in the CHB method, 
we could get the same values starting with the two initial spin configurations in contrast 



with the conventional MC method [^. Thus, we conclude that the CHB method, in fact, 
relieve the difficulty of the computer simulation of the ANNNI model. 

We have proposed a new update algorithm of the spin configuration and demonstrated its 
efficiency in the ANNNI model in two-dimensions. We should note again that the algorithm 
is a natural generalization of the conventional heat bath algorithm and applicable to various 
systems with short-range interactions. Since a large fluctuation of the spin configuration may 
occur for every MC sweep, the CHB method is particularly useful for studying equilibrium 
properties of complex systems such as spin-glasses . We also note that we may perform 
much more effective MC simulation by combining the CHB method with extended ensemble 
methods such as the exchange MC method 

The authors wish to thank Dr. T. Nakamura for valuable discussions. They also wish 
to give their thanks to Dr. S. Fujiki for showing them his unpublished data of the ANNNI 
model. This work was partly financed by a Grant-in-Aid for Scientific Research from the 
Ministry of Education, Science and Culture. The simulations were made partly on FACOM 
VPP500 at the Institute for Solid State Physics in University of Tokyo. 
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FIGURES 

FIG. 1. The MC sweep dependences of the order parameter M2 at T = 0.9Ji starting with two 

initial spin configurations of a paramagnetic phase (PARA) and the (2, 2) antiphase (AP) by (a) 
the conventional MC method and (b) the CHB method. 

FIG. 2. The temperature dependences of the average value of the order parameter (M2) starting 
with two initial spin configurations of a paramagnetic phase (PARA) and the (2, 2) antiphase (AP) 
by (a) the conventional MC method and (b) the CHB method. In the conventional MC method, 
the averages were taken over 2 x 10"^ — 4 x 10"^ and 1 x 10^ — 1.5 x 10^ MC sweeps for T > l.OJi 
and T < 0.9 Ji, respectively. On the other hand, in the CHB method, those were taken over much 
smaller numbers of MC sweeps, i.e., 5 x 10^ - 10 x 10^ and 1 x 10^ — 1.5 x 10^ MC sweeps for 
T > 1.0 Ji and T < 0.9 Ji, respectively. 
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FIG.2. F. Matsubara et al. 



